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Abstract 

The multidimensional gas-kinetic scheme for the Navier-Stokes equations under gravitational fields 
[J. Comput. Phys. 226 (2007) 2003-2027] is extended to resistive magnetic flows. The non- 
magnetic part of the magnetohydrodynamics equations is calculated by a BGK solver modified due 
to magnetic field. The magnetic part is treated by the flux splitting method based gas-kinetic theory 
[J. Comput. Phys. 153 (1999) 334-352 ], using a particle distribution function constructed in the 
BGK solver. To include Lorentz force effects into gas evolution stage is very important to improve 
the accuracy of the scheme. For some multidimensional problems, the deviations tangential to the 
cell interface from equilibrium distribution are essential to keep the scheme robust and accurate. 
Besides implementation of a TVD time discretization scheme, enhancing the dynamic dissipation 
a little bit is a simply and efficient way to stabilize the calculation. One-dimensional and two- 
dimensional shock waves tests are calculated to validate this new scheme. A three-dimensional 
turbulent magneto-convection simulation is used to show the applicability of current scheme to 
complicated astrophysical flows. 
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1. Introduction 

Plasma is a very common phase of matter in the universe, especially under astrophysical cir- 
cumstances. For instance, star can be regarded as very hot ball of plasma. The solar surface 
activities are thought to be tightly related to the interaction between turbulent convection and 
magnetic flux tubes. The accretion flows are always magnetized. Magnetic fields are also impor- 
tant in the interstellar medium of spiral, barred and irregular galaxies. These problems usually 
are complicated and high-nonlinear. Along with the fast growth of computer power, numerical 
computation is becoming a more and more popular way to study them. 

During the past decades, to meet the demands of industrial design and scientific computation, a 
lot of efforts have been made to develop accurate and robust numerical schemes for supersonic flows. 
Various methods, such as Godnov scheme, piecewise parabolic method (PPM) [3], total variation 



diminishing (TVD) scheme[6j] and gas-kinetic BGK scheme 2l|, |2__], have been invented based on 
different philosophies. When calculating the shock waves, we face the dilemma of keeping accuracy 
or robustness. In order to stabilize the calculation, we need dissipation. While an unnecessary 
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dissipation would reduce the accuracy of the scheme. Hence, the capability of a shock-capturing 
scheme depends on how smart the introduction of the numerical dissipation is. In the gas-kinetic 
BGK scheme, this is done by particle collisions. Due to collision, the particles will approach a 
thermally equilibrium state, e.g., Maxwellian distribution. Practice shows that this is a smart way 
to control dissipation. The resulted supersonic flow solver is accurate and robust. 

At the same time, some schemes for magnetohydro dynamics (MHD) have been developed or 
extended from their hydrodynamics (HD) version. For example, the PPM has been extended 
to mulitdimensional MHD by Dai and Woodward [4]. They also proposed an approximate MHD 
Riemann solver and an approach to maintain the divergence free condition of magnetic field [jj. 
A second-order accurate TVD Roe-type upwind MHD scheme and its multidimensional version 
were presented in 0, Most of them are based on characteristic decomposition. Because of the 
non-strictly hyperbolicity of the MHD system, it is hard to validate the MHD eigensystem. The gas- 



kinetic theory based flux splitting method for ideal MHD developed by Xu 19] is a robust and first- 
order scheme. It was extended to higher order and multidimensional by Tang and Xu 14|. Without 
wave decomposition, it is an efficient scheme. For the flux vector splitting (FVS) schemes [13], the 
flux function is split to F = F + + F~ , where the numerical dissipation is proportional to the 
numerical resolution. An equilibrium part due to particle collisions is introduced in fl9( ]. i.e., 
F = (1 — a)F e + a(F + + F~). a £ [0, 1] is an adjustable parameter used to control the particle 
collisions. In doing so, the unnecessary numerical viscosity is reduced. 

For the BGK-Navier-Stokes (BGK-NS) solver, the particle collisions are controlled by the so- 
lution of BGK equation[lt] instead of an arbitrary parameter. The reduction of dissipation only 
takes place where it is physically necessary. In oder to apply the BGK-NS solver to astrophysical 
flows, the magnetic field has to be included. The problem is that an implementation of the mag- 
netic terms with arbitrary descretization scheme might give rise to complicated truncation error 



behaviour. As analyzed by [17[, it is better to use the same difference operator to discretize all 
the terms in the equations, otherwise, the interaction between different kinds of truncation errors 
would introduce unphysical effect, e.g., numerical heating or numerical dispersion. In the current 
study, the non-magnetic part of MHD equations is solved by BGK-NS scheme, the magnetic part 
is solved by flux splitting method based on the solution of BGK equation. 

This paper is organized as follows. Section [2] describes the extension of BGK scheme to include 
magnetic field. In Section [3] numerical experiments are performed to validate the current scheme. 
Section H] discusses selected problems related to the current scheme and some general issues of 
developing MHD scheme. Section [5] is the conclusion. 

2. Numerical Scheme 

In this section, I describe a scheme for solving the following resistive MHD equations under 
gravitational field: 



dp/dt = -V-pv, (1) 

dpv/dt = -V- (pvv-BB)- Vptot + V-S + pf-V^), (2) 

dE/dt = -V ■ [(£ + ftot)? - ? • S + F d - • i) - B x i)(V x B)] + pw ■ (-V</)), (3) 

dB/dt = -V • (vB - Bv-r)VB) (4) 
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where, ptot = Pg + Pmi E = Ei + E m + Ef., p g is the gas pressure, p m = \B ■ B is the magnetic 
pressure, Ei is the internal energy, E m = ^B ■ B is the magnetic energy and E^ = \pv 2 is the 
kinetic energy. (— V</>) is the gravitational acceleration, Fa the diffusive heat flux, £ the viscous 
stress tensor, n the magnetic resistivity. All the other symbols have their standard meaning. The 
above equations can be written in a compact form: 
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where U = (p, pv,B,E) and Q = (0, p(— V</>), 0, pv ■ (— are cell averaged quantities. F x , F y 
and F z are defined on the cell interface. In each direction, the fluxes can be split into two parts. 
For example, F x can be further written as F x = F x b g k + i 7 ^, split; where 
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We discretize the non-magnetic part -F^^gk by gas-kinetic BGK scheme, magnetic part -F^sput by 
gas- kinetic theory based flux splitting method. 

2.1. A BGK-NS solver for the non-magnetic part 

The current scheme calculates i^^gk by the BGK-NS solver [13], which is needed to be modified 
to take into account the effects of magnetic field. First, the total pressure includes magnetic 
pressure, and total energy magnetic energy. Secondly, the effects of Lorentz force are added into 
the BGK equation[lj], i.e., 



df 



+ c • V/ + (-V0 + F Lorenzt ) ■ V s f 



9-f 



(8) 



dt t 
where f(t,x,y,z,c x ,Cy,c z ,£) is the particle distribution function defined on phase space, g is the 
Maxwellian distribution, 
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r is relaxation time, c= (c x ,c y ,c z ) is the particle micro velocity. £ is the internal variable, which 
has N internal degree of freedom. The value of N is determined by the configuration of molecules. 
A is related to the gas temperature by A = m/2kT, m is the molecules mass, k is the Boltzmann 
constant, and T is the temperature. FL orentz = B x (V x B)/p is the acceleration due to Lorentz 



force. Equation (|8|) has a local approximate solution at cell inter f ace [171]: 



/ = {(1 - e-* /r ) + (eT tlr (t + r) - r)[(a l n c n + b l n {-V<t> + F Lorentz ) n )H[c n ] 

+«c n + b r n (-V(/) + F Lorentz ) n )(l - H[c n }) 

+a t -ct + bf (-V0 + F Lorentz ) t ] + [t- r(l - e - t/T )]A}g 
+ e' l ' T { [l-af-ct-b 1 - (- V0 + F Lorentz )t 

-r(a l -c+P ■ (-V0 + FLorentz) + A l )]H[c n ]g l 
+ [1 - cf • ct - b r ■ (-V0 + F Lorentz )t 

-r((f • c + 6 r • (-V0 + FL orentz ) + ^ r )](l - H[c n ])g r }, (10) 

where the Heaviside function H[x] is defined by 

f 0,x < 0, 
[ 1, x > 0, 

x = is the location of cell interface. For a vector defined on the cell interface, its components 
are grouped into two categories, the component normal to the cell interface and those tangential 
components. For instance, for velocity c, c n is the normal component, are the tangential compo- 
nents. See Fig. [T]for an example, a" = [a l n ,a t ], a 1 = [a l n ,a l t ], b = [^,6 t ], V" = [6^,64], a = [a^,a[], 

W = [a^,a[], b = [bn,b t ], V = \b r n ,b r t ], A , A T and A are the coefficients measuring the deviations 
from Maxwellian due to various effects. Calculation of these coefficients is one of the key part of 
BGK scheme. The details are presented in [17j . The way that current scheme including magnetic 
field does not affect these calculations, go, g l , g^ are equilibrium and split Maxwellian, respectively. 
Once the solution of BGK is found, the fluxes of conservative quantities through the cell interface 
can be obtained, i.e., 



F x ,hgk = -J CxftpdS, (11) 

^ = [l,c,0,Q,0,\{S-c + e + B-B)} T , (12) 

where cG = dc x dc y dc z d^ is the volume element in phase space. The above integrals can be worked 
out analytically, which involves another key part of the BGK scheme - calculation of the moments 



of Maxwellian. The formula needed can be found in the appendix of [19( or [17j. 



2.2. A flux splitting solver for the magnetic part 

In order to perform the flux splitting, the solution of BGK equation obtained in previous section 
has to be split into left-side and right-side status, i.e., 



f(c n > 0) = {1 + (e~*/ r (t + r) - r)[(a l n c n + b l n (-V<t> + F Lorentz ) n ) 
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Figure 1: Two adjacent finite volumes - cells. (i,j) are the cell indices. Dashed lines are the cell interfaces. Fluxes of 
the conservative quantities, eg., F x and F y are calculated at the cell boundaries. The cell-averaged Uij is updated 
every time step. Solid lines show schematically the distribution of magnetic field constructed by van Leer limiter 
around cell interface if there is a shock. Arrows are a vector and its normal and tangential components. 

+a t -c t + b t - (-V0 + F Lorentz ) t ] + [t- r(l - e~ l l T )\A 
+ e-*/ T [a 1 -ct + l) ■ (-V0 + F Lorentz )t 

+r(a l -c + b 1 - (-V0 + F Lorentz ) + A l )]}g l , (13) 

f(c n < 0) = {1 + (e-^it + r) - r)[(a r n c n + b r n {-V<j> + F Lorentz ) n ) 
+a t -ct + b t - (-V0 + F Lorentz ) t ) + [t- r(l - e- t/r )}A 

+ e~ t/T [a r -ct + b r - (- V0 + F L orentz)t 

+T (^ . c - + . ( _ v</) + i?^^) + A r )M, (14) 

Two points need to be remarked on the above splitting. 

1. The tangential derivations due to derivatives and accelerations are very important for the 
BGK-NS solver. However, they do not affect the magnetic part evidently. This property 
seems case-dependent. For example, for the Orzang-Tang vortex problem, including tan- 
gential deviations is little bit dispersive. While for the cloud-shock interaction problem, it 
essentially improve the ability of the scheme. See Section [3] for the details. 

2. The right-hand side of Equation (fTUj) consists three parts, i.e., equilibrium part [• • -]go, split 
parts [• • -]g and [■ • -]gQ. The left and right parts introduce kinematic dissipation to stabilize 
calculation of discontinuity; while the equilibrium part reduces the kinematic dissipation 
to an appropriate value, thus improve the accuracy of the numerical scheme. These three 
terms together form a distribution which would approach an equilibrium state due to the 
particle collision. If we split the distribution (|10p into three parts accordingly, and replace 
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those squares with an adjustable constant a, (0 < a < 1), i.e., (1 — a)go, ag l and ag^, 
then the following fluxes construction is the same as in method 19)]. However, numerical 
experiment shows that this kind of splitting seems too dispersive, even without substituting 
those coefficients. This will be discussed further in Section [H 

Using the above split particle distribution, the magnetic part of fluxes can be split as: i^spiit = 
F l x + F x . F l x is computed as follows, 
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where (B l x , B l y , B l z ) and (B x , B y , B r z ) are the magnetic field at the left-hand side and right-hand 
side of the cell interface, respectively. They are reconstructed at the beginning of each time step 
at the cell interface. 

& = Bi( Xi ) + (x i+1/2 - Xi )Li, (16) 
B r = B i+1 (x i+ i) - (x i+ i - x i+ i/ 2 )L i+1 . (17) 
In the above expressions, the van Leer limiter is defined as, 

|s_i_||s_| 



Li — 5*(s+, S- 



+ s. 



where S(s + , s_) = sign(s+) +sign(s_), s + = (B i+1 - Bi)/(x i+l - Xi), and s_ = (B t - B i ^ 1 )/(x i - 
Xi-\). f} is constructed in the same way. Fig. [T] gives a schematic example for the reconstruction. 

F x is calculated in the similar way, using f and integrating on the interval [— oo, 0] . F y and 
F z are obtained by properly permuting indices. 

Some points need to be remarked on the above fluxes splitting, 

1. The advective magnetic terms are treated to be passively transported by the motions of 
particles. v y ,v z in ID problems and v z in 2D problems are also passively transported. Those 
non-advective terms are weighted by the particle distribution. 

2. Basically, the integrals on the right-hand-side of (|15p are separated according to the following 
rule. If a term explicitly contains the component of velocity norm to the cell interface, eg., 
v x , then v x is replaced by the particle velocity c x and grouped into the second integral. If 
not, it is grouped into the first integral. 
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3. An exception of the above rule is the term v l x B x B l x . Because in Eq. ([3D,we already include 
the magnetic pressure p m and magnetic energy density E m into the total pressure ptot and 
total energy E. The transported magnetic contributions to ptot and E, eg., v x B x B x should 
be canceled by the terms from BB ■ v in Eq. 0. The BGK scheme automatically divides p m 
and E m into two parts, eg., p m = (v x p m + c x p m )/2. In order to cancel these terms, we split 
v l x B l x B l x into: fc l x B l x B l x /2 and f l v l x B l x B l x /2, accordingly. 

2.3. Equation of State 

Internal variable £ is used to take account of the realistic effects of molecular structure to 
equation of state (EOS). For 3D diatomic molecule ideal gas, p£ 2 /2 is the sum of vibrational and 
rotational energy. £ is a function of the effective internal freedom degree K. For 3D cases, K = N . 
For 2D cases, K = N + 1, where the particle's motion in the additional third direction is included. 
For ID problems, K = N + 2. Instead of temperature T and internal energy Ei, in the gas 
evolution stage of the BGK scheme, A and internal freedom degree K are involved explicitly. For 
non-magnetic ideal gas, they are defined as: 



\ — P 

A " V 

K = A(Ei)X/p-3 
For general cases, we define effective A* and K* as following: 

A* - 9 
2ptot 

K* = A(E tot -E k )X*/p-3, 

where, ptot = Pg+Pr +Pm and E to t = E^ + E r + E m + E^. p m and E m are the pressure and energy 
due to magnetic field. p r and E r can be the pressure and energy due to radiation or ionization. 

2.4- Controlling Dissipation 

Besides MUSCL type of reconstruction, the numerical dissipation can be determined explicitly. 
The aforementioned collision time is defined as following, 

r = d At + C 2 ~ ^ ot | At + (22) 
\Ptot+Ptot\ Ptot 

where At is numerical time-step, \x is the dynamic viscosity, p\ ot , p\ ot are reconstructed total 
pressure at the left- and right-side of the cell interface. C2 = 1. For viscous flows, C\ = 0, for 
inviscid flows, C\ = 0.1 ~ 0.2. For magnetic flows, usually, the resistivity r\ is set to zero for 
ideal MHD. Sometimes, a small value of resistivity, e.g., rj = 0.008At is needed to stabilize the 
calculation. 

2.5. Ensuring V • B = 

The divergence free condition of magnetic field means the nonexistence of magnetic monopoles. 
Although in theory, they possibly exist, in reality they have never been found. On the contrary, 
caused by numerical truncation error, in MHD simulation, non-zero V • B is accumulated rapidly 
during time marching and destroys the physical consistence of the system. In practice, there are 
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(18) 
(19) 



(20) 
(21) 



different ways to suppress such kind of accumulation. A detailed comparison of these method is 
given by Toth [3]. In 8- wave scheme, the non-zero divergence is propagated away. More often, 
people use projection method to correct B obtained from MHD scheme at time-step n+1. Through 
B* = B n+1 + VG, where V 2 Q + V • B n+1 = 0, V • B* = is satisfied. Since this method needs 
solving a Possion equation, it is not very efficient for 3D problems solved on very large grids. In the 
current study, I use projection method for 2D programs. For 3D programs, I employ the so-called 
field- interpolated constrained transport scheme [5J to ensure the divergence free constraint. In this 
method, perfect cancellation of terms makes the divergence unchanging during time marching. It 
is simply and time saving. After obtaining the cell-averaged B at time i n+1 by gas kinetic scheme, 
we have magnetic field at time t n+1 ^ 2 , B n+1 / 2 = 0.5(i? n + B n+l ). B n+l / 2 is used to evaluate the 
right-hand-side of Eq. ([23]) . 

dB 

— = V x (v x B-rjV x B). (23) 
ot 

The left-hand-side of Eq. (|23p represents the values on the staggered mesh. After the values 
on staggered mesh are updated, the new cell-averaged values are obtained though interpolation. 
In principle, this method can ensure divergence free constraint to machine error. However, the 
truncation error occurring in boundary ghost cells propagates into the computational domain. 
Finally, we maintain V • B = to the truncation error. 

After the correction of the magnetic field, the total energy should be modified also to make the 
computation physically consistent, i.e., 

E n+1 = E n+1 _ I(£n+1 . jjn+l) + . ( 2 4) 



3. Numerical Tests 

For all the tests presented in this section, the van Leer limiter is used for the initial reconstruc- 
tion of the conservative variables. In order to make the comparison as precise as possible, for some 
problems, a constant numerical time-step is used. 

3.1. ID MHD Shock Waves 

To validate the scheme described in Section [2] and check the effects of Lorentz force in the gas 
evolution stage, I test the Biro-Wu ID MHD shock problem [2j. The calculations are performed 
on x G [0,2]. The initial conditions are set up with the left state (p,v x ,v y ,v z ,p g , B x , B y , B z ) = 
(1,0,0,0,1,0.75,1,0) and the right state (0.125,0,0,0,0.1,0.75,-1,0). Density p, gas pressure 
p g , total energy density E, x-velocity v x , y-velocity v y and y-magnetic field B y at time t = 0.2 
are plotted in Fig. [2] an d Fig. [3l When computing the solutions in Fig. [21 a second-order TVD 
Runge-Kutta scheme |1 01] is used to do time-marching. The solutions in Fig. [3]are obtained by Euler 
time-marching scheme. Fig. [2] compares the gas evolutions govern by the Maxwillian or the solution 
of BGK equation. We can see the stability is increased by using solution of BGK equation. Fig. [3] 
shows the effect of including the Lorentz force in gas evolution stage. The oscillations between the 
front of shock and the front of fast wave are greatly suppressed. This kind of suppressing is not 
obvious in second order time marching scheme, which means that the Lorentz force plays a role of 
first order time accuracy effects in the gas evolution. 
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Figure 2: The MHD shock tube test. A second order TVD Runge-Kutta scheme is used for the time stepping. 
Dots are the result of the Maxwellian based gas-kinetic flux splitting scheme with a = 0.9. Lines are the result of 
the current MHD BGK scheme. The simulations use 400 cells and a time step of At — 0.2Ax, corresponding to a 
Courant constant of ~ 0.78. A ratio of specific heats 7 = 2 is adopted. The plotted quantities are density p, gas 
pressure p g , total energy density E, x-velocity v x , y- velocity v y and y-magnetic field B y at t — 0.2. 
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Figure 3: The MHD shock tube test. The Euler time-stepping scheme is used. Dots are the result of the current 
MHD BGK scheme with effects of Lorentz force included in the gas evolution stage. Lines are the result of the MHD 
BGK scheme without considering Lorentz force in the particle distribution function. The simulations use 400 cells 
and a time step of At = 0.2Ax, corresponding to a Courant constant of 0.78. A ratio of specific heats 7 = 2 is 
adopted. The plotted quantities are density p, gas pressure p g , total energy density E, x-velocity v x , y-velocity v y 
and y-magnetic field B y at t = 0.2. 
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3.2. Spherical Explosion[22] 



This test is used to test the influence of strong magnetic field to the shock wave propagation. 
The computational domain is [0, 1] x [0, 1]. Initially, the density is 1 everywhere. There is a high 
pressure region around the center with a radius of 0.1. The pressure inside and outside the central 
region are 100 and 1, respectively. The magnetic field is initialized as (B x , B y , B z ) = (0, 50/^^", 0). 
The ratio of specific heats 7 equals 2. A uniform 100 x 100 mesh is used for this problem. The 
boundaries are periodic. 

Fig. [3] shows a snapshot of the solution at t = 0.0105. Density, gas pressure, magnetic pressure 
and kinetic energy are visualized with 25 contours for each quantities. Due to the strong magnetic 
field, the spherical explosion are highly anisotropic. The shocks also form in the magnetic field. 
Fig. compares the density distributions from different scheme. These distributions are taken 
along y = 0.5 at time t = 0.0.105. Solid line is obtained by the current BGK MHD scheme. Dots 
are from Maxwellian based gas-kinetic scheme with a = 0.95. The post-shock region indicates that 
the current scheme is less dissipative than the Maxwellian based scheme. 

3.3. Orszag-Tang Turbulent Vortex 

The 2D Orszag-Tang turbulent vortex problem [3] is widely used by computational fluid dy- 
namics (CFD) community to test the MHD scheme because of its complicated interaction between 
different waves generated as the vortex system evolving. The computations are preformed on a 
domain of [0 : 2tt] x [0 : 2ir\. The initial conditions are 

p(x,y) =>y 2 ,u x = - sin (y),u y = sin(x), 
p g (x,y) = -y,B x = - sin (y), By = sin(2x). (25) 

Fig. [6] shows a snapshot of the solution at time t = 3. Twenty five contours for density, gas 
pressure, magnetic pressure, kinetic energy are plotted. This solution is obtained by BGK MHD 
scheme with a second-order TVD Runge-Kutta time marching scheme. The effects of tangential 
derivatives are excluded during gas evolution. When considering the acceleration due to Lorentz 
force, only the normal component is included in the particle distribution function. The divergence 
free condition of magnetic field is ensured by projection method. Around this time, the four spiral 
arms like shock fronts are propagating anticlockwise. 

Fig. E] compares the density distributions from different splitting methods. A solution of 
Maxwellian based gas-kinetic splitting scheme is also displayed in panel (a). In (b), the Lorzent 
force is switched off in gas evolution stage. In (c), the tangential deviations are switched on. In 
(d), the tangential component of Lorentz force is included. The oscillations caused by tangential 
deviations are evident. The effect of Lorentz force is not obvious in (b). 

A more detailed comparison is given in Fig. [8j in which the gas pressure distributions along 
y = 0.6257T at t = 3 are plotted. At this time, these waves are propagating forward along x 
axis. There are three obvious differences in these distributions. The most obvious different is the 
altitude of the hump post first shock - the shock at x ~ 4.3. This hump goes through a long 
term evolving. This difference is caused by particle distributions. Fig. [5J indicates that the BGK 
solution based scheme is less dissipative than the Maxwellian based scheme. The difference in the 
fluctuations behind second shock (the one locates at x ~ 1.6) is mainly caused by including the 
Lorentz force effects in the gas evolution. The fluctuations obtained by current scheme is more 
smooth. The third apparent difference occurring at x ~ 5.25 is caused by the implementation of 
tangential deviations in particle distribution functions. 
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Figure 4: A snapshot of the solution to the spherical explosion problem at time t = 0.0105. The calculation is 
performed on a uniform mesh of 100 x 100 grids, using the current BGK MHD scheme. 25 contours for each 
quantity are plotted. Numerical time step is fixed at At = 0.008, corresponding to a Courant number ~ 0.67. 
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Figure 5: The density distributions of the spherical explosion problem along y = 0.5 at time t = 0.0105. Solid line is 
obtained by the current BGK MHD scheme. Dots are from the Maxwellian based gas-kinetic scheme with a — 0.95. 

3-4- Cloud-Shock Interaction^] 

This simplified astrophysical problem describes the disruption of a denser cloud by a magne- 
tosonic shock. The computational domain is 1 x 1 box solved on a 200 x 200 mesh. Initially, 
the rightward-propagating shock locates at x = 0.5 with a Mach number 10. The left and right 
states are (p, v x , v y , v z ,p g , B x , B y , B z ) = (3.86859,11.2536,0,0,167.345,0,2.1826182,-2.1826182) 
and (1, 0, 0, 0, 1, 0, 0.56418985, 0.56418985), respectively The circular cloud is centered at (0.25, 0.5) 
with a radius 0.15 and density 10. It is in hydrostatic with its surrounding. 

This problem is solved by our 3D BGK MHD code with a very thin third dimension. The 
divergence free condition is ensured by the field-interpolated constrained transport scheme. The 
Euler scheme is used for time-stepping. If the magnetic resistivity is set rj = 0, the calculation easily 
crashes. Instead of employing the 2nd order TVD Runge-Kutta scheme and the time-consuming 
projection method ensuing V • B = 0, I increase the resistivity slightly. Numerical experiments 
show a r\ = 0.008A£ and C\ = 0.2 are enough to keep the computation stable. The tangential 
deviations in the distribution function plays a very important role for this case. If we switch 
off those tangential deviations, the oscillations occur widely and crash the computation. Fig. [10] 
shows the distributions of various quantities along the middle line y = 0.5. It is clear that the 
current scheme can capture shock with 4 ~ 7 cells, even with the artificially enhanced viscosity 
and resistivity. 

3.5. 3D Turbulent Magneto-convection 

The 3D turbulent magneto-convection in the stellar interior is calculated to show the appli- 
cability of the current BGK MHD scheme to complicated astrophysical flows. This is one of the 
most difficult MHD problem in astrophysics. It is essentially a multi- length-scale and time-scale 
problem, which makes its numerical modeling very time- and memory- consuming. In the deep 
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Figure 6: A snapshot of the solution to the Orszag-Tang MHD turbulent problem at time t — 3. The calculation 
is performed on a uniform mesh of 192 x 192 grids, using BGK MHD scheme without including the tangential 
derivatives in the gas evolution stage. 25 contours for each quantity are plotted. Numerical time step is fixed at 
At = 0.008, corresponding to a Courant number « 0.67. 
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Figure 7: The density distributions of the Orszag-Tang problem at time t — 3. (a) Maxwellian based gas-kinetic flux 
splitting with a = 0.7; (b) current scheme with Lorentz force switched off in the gas evolution stage; (c) tangential 
derivates switched on; (d) tangential component of Lorentz force switched on. 
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Figure 8: Gas pressure distributions along y — 0.625-7T at t — 3. Dots are from the Maxwellian based gas kinetic flux 
splitting scheme. Lines are from BGK MHD scheme, (a)tangential deviations switched off; (b)Lorentz force switched 
off in the gas evolution; (c)tangential derivatives switched on; (d) tangential component of Lorentz force switched 
on. 
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Figure 9: The solution of cloud-shock interaction problem at time t = 0.055. The solution is obtained by 3D BGK 
MHD code with a very thin third dimension. Density, gas pressure magnetic pressure and kinetic energy are plotted. 
V • B = is ensured by the field-interpolated constrained transport method. A 200 x 200 mesh and Curant number 
0.7 are used. 
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Figure 10: The distributions of various quantities along y = 0.5 for cloud-shock interaction problem at time t = 0.055. 
In the left panel: dots are tot pressure; solid line is gas pressure. In the right panel: dots are magnetic pressure; 
solid line is density. 



stellar atmosphere, the convective motions are subsonic. However just beneath the upper radiative 
zone, there is thin highly superadiabatic layer, where complex waves, including supersonic waves 
generated. The computation is easily interrupted if the code is not robust, especially, during the 
thermal relaxation phase of the model. 

Here I calculate a 3D sandwich model of stellar convection: one convective layer locates in 
between two stable layers. The stable layers are radiative. Radiation transfer is treated by diffusion 
approximation. Under constant gravity, the system is initially in hydrostatic state and vertically 
stratified. The method to construct the initial stratification can be find in The middle 

layer is slightly superadiabatic. The stable layers are subadiabatic. The vertical depth extents 
5.3 pressure scale heights. The sub-grid scale turbulent behaviors are mimicked by Smagorinsky 



model Il6l |. The aspect ratio is 3/1 (width/depth). The horizontal boundaries are periodic. 



The upper boundary and lower boundary are impenetrable. All quantities are scaled such that the 
density, pressure, temperature at the top boundary are unit. The length is scaled by the depth 
of the computation domain. This initial state will undergo an adjustment, significantly near the 
vertical boundaries and stable-unstable interface, moderately in the convective region. Finally, the 
system will be thermally relaxed and approach an dynamic equilibrium state, indicated by the 
balance between outgoing energy flux through top and input energy flux from bottom. 

When the non-magneto convection nearly approaches the relaxed state, I impose a uniform 
vertical magnetic field to a snapshot of the numerical solutions. Its strength equals the equipartition 
value near the top. At the top and bottom, the vertical field boundary condition is adopted for 
the magnetic field, i.e., 

B x = B y = 0, ^ = 0. (26) 
The computation was run on the supercomputer with a 336 threading code parallelized by massage 
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Figure 11: Results of 3D turbulent magneto-convection. Horizontal slices at z = 0.82. Upper panel: color contours 
are the vertical component of the velocity; arrows are the horizontal components. Lower panel: color contours are 
the magnetic strength; arrows are the horizontal components of the magnetic field. 
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Figure 12: Results of 3D turbulent magneto-convection. Vertical slices at y = 1.67. Upper panel: color contours 
are the temperature fluctuation from horizontal mean; arrows are the projection of the velocity on x-z plane. Lower 
panel: color contours are the magnetic strength; arrows are the projection of the magnetic field on x-z plane. 
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Figure 13: Results of 3D turbulent magneto-convection. Time history of maximum Mach number, t — is the 
initialization time of magentic field. Solid line indicates the property of non-magnetic convection. Dotted line 
indicates the magnetic convection 



passing interface (MPI) protocol. The computation of non-magnetic convection took ~ 36 days, it 
is still not completely relaxed. The upper boundary can transport around 85 percent of the input 
energy fluxes. The magneto convection was run for ~ 7 days. At the end, 93 percent input fluxes 
can be transported out from upper lid. 

Fig. [TT] shows the projections of 3D velocity and magnetic field on a horizontal x-y plane, at 
z = 0.82. Upper panel shows the velocity field. Lower panel shows the magnetic field. Color 
contours are the vertical component. Arrows are horizontal fields. Fig. [T2] displays vertical slices 
of the flows. In upper panel, the 2D velocity field is imposed on the fluctuation of temperature 
from the horizontal mean. In lower panel, the 2D magnetic field is imposed on the contour of 
magnetic strength. Fig. [13] shows the time history of the maximum Mach number of the system. 
The magnetic field suppresses the motion evidently. Before superposing magnetic field, the highest 
speed motion is nearly supersonic. 

4. Discussion 

The current scheme can be regarded as a hybrid method, a combination of BGK-NS solver and 
flux splitting scheme. Without introducing parameter-controlled equilibrium state, the treatment 
of the magnetic part of fluxes is similar to the FVS scheme. Unlike the FVS scheme and the gas- 
kinetic MHD method developed by Xu[19], where the Maxwellian distribution is used, the current 
method splits the solution of BGK equation. The collision mechanism is already included in the 
split distribution. Physically, the current scheme is supposed to be better. However, numerical 
tests show that the improvement is not very great. For some cases, it is dispersive. 
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For non-magnetic flows, the tangential deviations from equilibrium state in particle distribution 
function play a role in improving the accuracy of the scheme. After implementing the magnetic 
field, the effects become case-dependent. The effects of tangential terms can be thought as reducing 
the numerical dissipation. For some cases, eg., 2D Orzang-Tang vortex, this kind of reduction seems 
too much, even with a 2nd order TVD Runge-Kutta time marching scheme. To get around with this 
kind of situation, it is better to switch off these tangential deviations terms in particle distribution 
functions. 

As stated before, the particle solution from solving BGK equation can be split into three parts, 
i.e., [■ • -]go, split parts [• • -]g l and [• • -J^g, just like Maxwellian based scheme, i.e., (1 — a)go, ag l and 
ag^. The adjustable parameter is replaced by complicated coefficients from solving BGK equation. 
It sounds more physical. Numerical tests again show that this kind of consideration would make 
a very dispersive scheme. This is because in Maxwellian based scheme, the particle collisions are 
controlled by an universal constant a, so the numerical dissipation from MUSCL-type reconstruc- 
tion is added in whole computational domain. In the BGK scheme, the numerical dissipation 
only become large near the region where the particle distribution deviates from equilibrium state 
significantly. More specifically, the equilibrium term and non-equilibrium terms from (|10p can be 
expressed as [1 + e~*/ T + • • -]go, split parts [e~*/ T (- • -)]g l and [e -i / T (- • -)]<7o- e - */ 7 " is a very small 
number, because usually we have r < t/10, which means the split parts can be ignored, r only 
becomes comparable to t near the high non-equilibrium region. This is reason why BGK is a 
very robust and accurate scheme for non-magnetic shock wave problems. For magnetic flows, the 
equilibrium part [1 + e~*/ T + • • -]go needs to be split also to include enough dissipation to smooth 
the dispersion. 

Fig. [3] shows clearly that the dispersion occurs near the front of fast wave and propagates 
backward. Microscopically, transport of flow properties can be regarded as a consequence of particle 
motions and the amount of flux is proportional to the number of involved particles. This is not 
true for the magnetic field, even those advective terms. This is because the unlike charge or mass, 
there is no specific amount magnetic is carried by each particle. If we treat the non-magnetic part 
by a second order time scheme, while the magnetic field by a first oder time scheme. Although, 
the magnetic terms are weighted by second order particle distributions, the difference of time 
accuracy of these two parts is still big. The truncation errors propagate at different speed, which 
amplifies the numerical dispersion. A second-order TVD Runge-Kutta scheme can smear out these 
oscillations. 

The second-order TVD time scheme usually catches a sloping shock front. A way to improve this 
is to design a very smart parameter a that can introduce numerical dissipation both near the shock 



and the fast waves. In 15[ a is adaptive according to the discontinuity in magnetic pressure. A 



more consistent way to construct a gas-kinetic theory based MHD scheme is to considers separately 
the distribution of charged particle (electrons and ions) and electrically neutral particle. Model 
the long distance interactions between charged particles and then integrate the electric field and 
magnetic properties from the motion of these charged particles. In such case, the EOS play a 
very important role, since it determines the number of charged particles. Obviously, this is a very 
complicated method. A more detailed discussion is beyond the scope of current study. 

The gas-kinetic MHD scheme developed in is an efficient method. For the 2D Orzang- 
Tang vortex problem tested in this paper, it can be two times faster than the current scheme 
(with ~ 0.47 time-consuming). However, in application to astrophysical problems, a lot of things, 
such as radiation transfer, turbulent modeling, moving mesh generating, need to be implemented. 
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The hydrodynamics part of the code is not the major agent consuming the computation time. 
Sometimes, even the ensuring divergence free constraint takes a lot time. For instance, for the 2D 
tests in this paper, a Possion equation is solved to ensue V ■ B = 0. This involves an implicit solver. 
The time-consuming of an implicit solver is generally not linearly proportional to the grid points. 
The MHD code based on current scheme is capable of calculating the complicated HD problem 
without extra efforts. At the same time, it is also a high order MHD solver. 

5. Conclusion 

This paper presents a feasible way to extend the multidimensional gas-kinetic BGK scheme to 
MHD problems. The non-magnetic part is solved by the BGK scheme modified due to magnetic 
field. While the magnetic part is calculated by gas-kinetic flux splitting based on a solution of BGK 
equations. Besides keeping high robustness and accuracy for HD problem, the current scheme is 
also a high order MHD solver. The effects of tangential deviations and Lorentz force in particle 
distribution are numerically studied. Although case-dependent, it is better to include them for 
most of the cases. Numerical tests also show that compared to the Maxwellian based gas- kinetic 
flux splitting scheme, the current scheme is more stable for ID problems and less dissipative for 
multidimensional problems. The 3D turbulent magneto-convection test indicates the applicability 
of current scheme to complicated astrophysical flows. The gas- kinetic theory based flux splitting 
methods for MHD problems can be improved by designing a more smart parameter a. A full 
gas-kinetic scheme for MHD needs to consider the complicated physical processes, such as particle 
charging and long distance Coulomb force. 
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